#include "StdAfx.h"
#include "MultiscaleNet.h"
#include <iterator>
#include <cmath>
#define edistance(x1, x2, y1, y2) ((double)x1-(double)x2)*((double)x1-(double)x2)+((double)y1-(double)y2)*((double)y1-(double)y2)
#define naturalnumber 2.7182818284590452353602874713527

//Constructor
MultiscaleNet::MultiscaleNet(unsigned int sizex, unsigned sizey, unsigned int sizez, unsigned int numberofdata, unsigned maxcycle,
							 std::vector<double> maxvalue, std::vector<double> minvalue)
{
	//Assign the attributes of the SOM
	srand(static_cast<unsigned int>(time(NULL)));
	sizeofsom.x = sizex;
	sizeofsom.y = sizey;
	sizeofsom.z = sizez;
	datadimension = numberofdata;
	trainingcycle = maxcycle;
	mapradius = min(sizex, sizey)/4;
	initialrate = (double)1.0;
	finalrate = (double)0.1;
	initialsize = (double)mapradius;
	finalsize = (double)1.0;
	learningratescale = log(initialrate/finalrate);
	neighborhoodscale = log(initialsize/finalsize);
	currentiteration = 0;
	currentcycle = 0;

	//Create the map
	somoutputlayer.resize(sizeofsom.x);
	for(unsigned int i=0; i<sizeofsom.x; i++)
	{
		somoutputlayer[i].resize(sizeofsom.y);
		for(unsigned int j=0; j<sizeofsom.y; j++)
		{
			somoutputlayer[i][j].resize(sizeofsom.z);
			for(unsigned int k=0; k<sizeofsom.z; k++)
			{
				_3dpoint currentposition = {i, j, k};
				somoutputlayer[i][j][k].position = currentposition;
				somoutputlayer[i][j][k].weight = MultiscaleNet::GenerateDoubleVector(maxvalue, minvalue);
				if(somoutputlayer[i][j][k].weight.size()==0)
					throw string("Size of max or min not correct");
			}
		}
	}

	return;
}

//Create a multiscale net without defining the max or min weights
MultiscaleNet::MultiscaleNet(unsigned int sizex, unsigned int sizey, unsigned int sizez, unsigned int numberofdata, unsigned int maxcycle)
{
	//Assign the attributes of the SOM
	srand(static_cast<unsigned int>(time(NULL)));
	sizeofsom.x = sizex;
	sizeofsom.y = sizey;
	sizeofsom.z = sizez;
	datadimension = numberofdata;
	trainingcycle = maxcycle;
	mapradius = min(sizex, sizey)/4;
	initialrate = (double)1.0;
	finalrate = (double)0.1;
	initialsize = (double)mapradius;
	finalsize = (double)1.0;
	learningratescale = log(initialrate/finalrate);
	neighborhoodscale = log(initialsize/finalsize);
	currentiteration = 0;
	currentcycle = 0;

	//Create the map
	somoutputlayer.resize(sizeofsom.x);
	for(unsigned int i=0; i<sizeofsom.x; i++)
	{
		somoutputlayer[i].resize(sizeofsom.y);
		for(unsigned int j=0; j<sizeofsom.y; j++)
		{
			somoutputlayer[i][j].resize(sizeofsom.z);
			for(unsigned int k=0; k<sizeofsom.z; k++)
			{
				_3dpoint currentposition = {i, j, k};
				somoutputlayer[i][j][k].position = currentposition;
				somoutputlayer[i][j][k].weight = MultiscaleNet::GenerateDoubleVector(datadimension);
				if(somoutputlayer[i][j][k].weight.size()==0)
					throw string("Size of max or min not correct");
			}
		}
	}

	return;
}

//Create the multiscale net from a file
MultiscaleNet::MultiscaleNet(std::string filepath)
{
	try
	{
		MultiscaleNet::LoadNetFromFile(filepath);
	}
	catch(string exceptionstring)
	{
		throw exceptionstring;
	}
}

//Generate a random double vector
vector<double> MultiscaleNet::GenerateDoubleVector(std::vector<double> maxvalue, std::vector<double> minvalue)
{
	vector<double> doublevector;
	//If the max and min size are not correct, return an empty double vector
	if(maxvalue.size()!=minvalue.size()&&maxvalue.size()!=sizeofsom.z)
	{
		return doublevector;
	}

	doublevector.resize(maxvalue.size());
	for(unsigned int i=0; i<doublevector.size(); i++)
	{
		//Generate random doubles between 0~1
		doublevector[i] = (double)(abs(rand())+1.0)/(double)(abs(rand())+1.0);
		if(doublevector[i]>1.0)
			doublevector[i]=(double)1.0/doublevector[i];
		//Rescale the random doubles to minvalue~maxvalue
		doublevector[i] = doublevector[i]*maxvalue[i]+((double)1.0-doublevector[i])*minvalue[i];
	}
	return doublevector;
}

//Generate a random double vector
vector<double> MultiscaleNet::GenerateDoubleVector(unsigned int numberofdata)
{
	vector<double> doublevector;

	doublevector.resize(numberofdata);
	for(unsigned int i=0; i<doublevector.size(); i++)
	{
		doublevector[i] = (double)(abs(rand())+1.0)/(double)(abs(rand())+1.0);
	}
	return doublevector;
}

//Compute the distance between two vectors
double MultiscaleNet::Distance(std::vector<double> vector1, std::vector<double> vector2)
{
	//Validate the size
	if(vector1.size()!=vector2.size())
		throw string("Size of vector for distance measure not equal");
	unsigned int vectorsize = (unsigned int)vector1.size();

	//Compute the distance
	double distance = 0;
	for(unsigned int i=0; i<vectorsize; i++)
		distance += pow(vector1[i]-vector2[i], (double)2.0);

	return distance;
}

//Find the winners
_3dpoint MultiscaleNet::FindWinner(vector<vector<double>> input)
{
	double mindistance = (double)ULONG_MAX;
	_3dpoint winnerposition = {UINT_MAX, UINT_MAX, UINT_MAX};

	if(input.size()>sizeofsom.z)
		throw string("Scale of input large than layers of SOM");

	for(unsigned int i=0; i<sizeofsom.x; i++)
	{
		for(unsigned int j=0; j<sizeofsom.y; j++)
		{
			unsigned int numberofloop = sizeofsom.z - (unsigned int)input.size();
			for(unsigned int k=0; k<numberofloop; k++) //Start matching from 0 until numberofloop
			{
				double currentdistance = 0;
				for(unsigned int l=0; l<input.size(); l++)
				{
					//Sum the distance at each scale up
					try
					{
						currentdistance += Distance(input[l], somoutputlayer[i][j][l+k].weight)/pow(2.0, (double)l);	/*Different weights can be placed here in order to consider the effects at different scale*/
					}
					catch(string exception)
					{
						throw exception;
					}
					//Compare with the current min distance
					if(currentdistance<mindistance)
					{
						mindistance = currentdistance;
						winnerposition = somoutputlayer[i][j][k].position;
					}
				}
			}
		}
	}

	//If no winner is found
	if(winnerposition.x == UINT_MAX)
		throw string("No winner found");

	return winnerposition;
}

//Find the winners and storing the distance
_3dpoint MultiscaleNet::FindWinner(std::vector<vector<double> > input, double &distance)
{
	double mindistance = (double)ULONG_MAX;
	_3dpoint winnerposition = {UINT_MAX, UINT_MAX, UINT_MAX};

	if(input.size()>sizeofsom.z)
		throw string("Scale of input large than layers of SOM");

	for(unsigned int i=0; i<sizeofsom.x; i++)
	{
		for(unsigned int j=0; j<sizeofsom.y; j++)
		{
			unsigned int numberofloop = sizeofsom.z - (unsigned int)input.size();
			for(unsigned int k=0; k<numberofloop; k++) //Start matching from 0 until numberofloop
			{
				double currentdistance = 0;
				for(unsigned int l=0; l<input.size(); l++)
				{
					//Sum the distance at each scale up
					try
					{
						currentdistance += Distance(input[l], somoutputlayer[i][j][l+k].weight)/pow(2.0, (double)l);	/*Different weights can be placed here in order to consider the effects at different scale*/
					}
					catch(string exception)
					{
						throw exception;
					}
					//Compare with the current min distance
					if(currentdistance<mindistance)
					{
						mindistance = currentdistance;
						winnerposition = somoutputlayer[i][j][k].position;
					}
				}
			}
		}
	}

	//If no winner is found
	if(winnerposition.x == UINT_MAX)
		throw string("No winner found");

	distance = mindistance;
	return winnerposition;
}

//Update the weights
void MultiscaleNet::UpdateWeights(_3dpoint winner, vector<vector<double>> input)
{
	unsigned int startx = (unsigned int)max((int)winner.x-(int)neighborhoodsize, (int)0);
	unsigned int starty = (unsigned int)max((int)winner.y-(int)neighborhoodsize, (int)0);
	unsigned int endx = min(sizeofsom.x-1, winner.x+(unsigned int)neighborhoodsize);
	unsigned int endy = min(sizeofsom.y-1, winner.y+(unsigned int)neighborhoodsize);
	unsigned int startz = winner.z;
	unsigned int endz = min((unsigned int)input.size()+winner.z-1, sizeofsom.z-1);
	//Update the weights at each hierarchy
	for(unsigned int i=startz; i<=endz; i++)
	{
		for(unsigned int j=startx; j<=endx; j++)
		{
			for(unsigned int k=starty; k<=endy; k++)
			{
				double distancetowinner = edistance(winner.x, j, winner.y, k);
				//If the node is inside the neighborhood function
				if(distancetowinner<=(neighborhoodsize*neighborhoodsize))
				{
					if(input[i-startz].size()!=somoutputlayer[j][k][i].weight.size())
						throw string("Size of node weight and input not equal");
					//Compute the neighhood function and update the weights
					double neighborhood = pow(naturalnumber, -distancetowinner/(neighborhoodsize*neighborhoodsize*(double)2.0));
					unsigned int sizeofnodeweight = (unsigned int)input[i-startz].size();
					for(unsigned int l=0; l<sizeofnodeweight; l++)
						somoutputlayer[j][k][i].weight[l] += learningrate * neighborhood * (input[i-startz][l]-somoutputlayer[j][k][i].weight[l]);
				}
			}
		}
	}

	return;
}

//update the weights including the patterns
void MultiscaleNet::UpdateWeights(_3dpoint winner, std::vector<vector<double> > input, unsigned int inputpattern)
{
	unsigned int startx = (unsigned int)max((int)winner.x-(int)neighborhoodsize, (int)0);
	unsigned int starty = (unsigned int)max((int)winner.y-(int)neighborhoodsize, (int)0);
	unsigned int endx = min(sizeofsom.x-1, winner.x+(unsigned int)neighborhoodsize);
	unsigned int endy = min(sizeofsom.y-1, winner.y+(unsigned int)neighborhoodsize);
	unsigned int startz = winner.z;
	unsigned int endz = min((unsigned int)input.size()+winner.z-1, sizeofsom.z-1);
	//Update the weights at each hierarchy
	for(unsigned int i=startz; i<=endz; i++)
	{
		for(unsigned int j=startx; j<=endx; j++)
		{
			for(unsigned int k=starty; k<=endy; k++)
			{
				double distancetowinner = edistance(winner.x, j, winner.y, k);
				//If the node is inside the neighborhood function
				if(distancetowinner<=(neighborhoodsize*neighborhoodsize))
				{
					if(input[i-startz].size()!=somoutputlayer[j][k][i].weight.size())
						throw string("Size of node weight and input not equal");
					//Compute the neighhood function and update the weights
					double neighborhood = pow(naturalnumber, -distancetowinner/(neighborhoodsize*neighborhoodsize*(double)2.0));
					unsigned int sizeofnodeweight = (unsigned int)input[i-startz].size();
					for(unsigned int l=0; l<sizeofnodeweight; l++)
						somoutputlayer[j][k][i].weight[l] += learningrate * neighborhood * (input[i-startz][l]-somoutputlayer[j][k][i].weight[l]);
					outputpatterncount[j][k][inputpattern] ++;
				}
			}
		}
	}

	return;
}

//Train the SOM
void MultiscaleNet::Training(vector<vector<vector<double>>> trainingdata, vector<_3dpoint>& trainingwinners,
							 vector<vector<vector<double>>>& trainingweights)
{
	//Resize the training results
	trainingwinners.resize(trainingdata.size());
	trainingweights.resize(trainingdata.size());
	unsigned int numberofdata = (unsigned int)trainingdata.size();
	unsigned int totaliteration = numberofdata*trainingcycle;
	trainingdistance.resize(trainingdata.size());

	//Start the training
	for(unsigned int i=0; i<numberofdata; i++)
	{
		//Update the parameters
		learningrate = initialrate * pow(naturalnumber, -(double)(currentiteration)/(double)totaliteration*learningratescale);
		neighborhoodsize = initialsize * pow(naturalnumber, -(double)(currentiteration)/(double)totaliteration*neighborhoodscale);
		currentiteration ++;

		//Varify the current cycle
		currentcycle = currentiteration / numberofdata;
		if(currentcycle>trainingcycle)
			throw string("Max training cycle exceeded");

		//Find the winner and update the weights
		try
		{
			trainingwinners[i] = FindWinner(trainingdata[i], trainingdistance[i]);
			UpdateWeights(trainingwinners[i], trainingdata[i]);
		}
		catch(string message) //Throw the error message
		{
			throw message;
		}
		unsigned int currentnumberofscale = (unsigned int)trainingdata[i].size();
		trainingweights[i].resize(sizeofsom.z);
		for(unsigned int j=0; j<sizeofsom.z; j++)
			trainingweights[i][j] = somoutputlayer[trainingwinners[i].x][trainingwinners[i].y][trainingwinners[i].z].weight;
	}

	return;
}

//Train the SOM with known patterns
void MultiscaleNet::Training(std::vector<vector<vector<double> > > trainingdata, std::vector<unsigned int> inputpattern, 
							 std::vector<_3dpoint> &trainingwinners, std::vector<vector<vector<double> > > &trainingweights)
{
	//Resize the training results
	trainingwinners.resize(trainingdata.size());
	trainingweights.resize(trainingdata.size());
	unsigned int numberofdata = (unsigned int)trainingdata.size();
	unsigned int totaliteration = numberofdata*trainingcycle;
	trainingdistance.resize(trainingdata.size());

	//get the number of patterns in the set
	unsigned int numberofpatterns = *max_element(inputpattern.begin(), inputpattern.end());
	outputpatterncount.resize(sizeofsom.x);
	for(vector<vector<vector<int> > >::iterator i=outputpatterncount.begin(); i!=outputpatterncount.end(); ++i)
	{
		(*i).resize(sizeofsom.y);
		for(vector<vector<int> >::iterator j=(*i).begin(); j!=(*i).end(); ++j)
			(*j).insert((*j).end(), ((numberofpatterns+1)-(*j).size()), -1);
	}

	//Start the training
	for(unsigned int i=0; i<numberofdata; i++)
	{
		//Update the parameters
		learningrate = initialrate * pow(naturalnumber, -(double)(currentiteration)/(double)totaliteration*learningratescale);
		neighborhoodsize = initialsize * pow(naturalnumber, -(double)(currentiteration)/(double)totaliteration*neighborhoodscale);
		currentiteration ++;

		//Varify the current cycle
		currentcycle = currentiteration / numberofdata;
		if(currentcycle>trainingcycle)
			throw string("Max training cycle exceeded");

		//Find the winner and update the weights
		try
		{
			trainingwinners[i] = FindWinner(trainingdata[i], trainingdistance[i]);
			UpdateWeights(trainingwinners[i], trainingdata[i], inputpattern[i]);
			//UpdateWeights(trainingwinners[i], trainingdata[i]);
			//outputpatterncount[trainingwinners[i].x][trainingwinners[i].y][inputpattern[i]]++;
		}
		catch(string message) //Throw the error message
		{
			throw message;
		}
		unsigned int currentnumberofscale = (unsigned int)trainingdata[i].size();
		trainingweights[i].resize(sizeofsom.z);
		for(unsigned int j=0; j<sizeofsom.z; j++)
			trainingweights[i][j] = somoutputlayer[trainingwinners[i].x][trainingwinners[i].y][trainingwinners[i].z].weight;
	}

	//determine the patterns for each node
	if(!outputpatterncount.empty())
	{
		outputpattern.resize(outputpatterncount.size());
		vector<vector<unsigned int> >::iterator outputite = outputpattern.begin();
		for(vector<vector<vector<int> > >::iterator i=outputpatterncount.begin(); i!=outputpatterncount.end(); ++i)
		{
			(*outputite).reserve((*i).size());
			for(vector<vector<int> >::iterator j=(*i).begin(); j!=(*i).end(); ++j)
			{
				vector<int>::iterator maxiterator = max_element((*j).begin(), (*j).end());
				if(*maxiterator == -1)		//if the max ite is -1 (no winner has been mapped to the node)
					(*outputite).push_back(-1);			//push -1 pattern into the node
				else
					(*outputite).push_back(static_cast<int>(distance((*j).begin(), maxiterator)));	//push the pattern into the node
			}
			++outputite;
		}
	}

	return;
}

//Get testing results
void MultiscaleNet::GetTestWinnerWeight(vector<vector<double>> testinput, _3dpoint& position, vector<vector<double>>& weight)
{
	try
	{
		position = FindWinner(testinput);
	}
	catch(string message)	//If winners not found
	{
		throw message;
	}

	//Get the winner and weights
	weight.resize(testinput.size());
	for(unsigned int i=0; i<weight.size(); i++)
		weight[i] = somoutputlayer[position.x][position.y][position.z+i].weight;
	
	return;
}

//Get testing results
void MultiscaleNet::GetTestWinnerWeight(vector<vector<double> > testinput, _3dpoint& position, vector<vector<double> >& weight, unsigned int& pattern)
{
	try
	{
		position = FindWinner(testinput);
	}
	catch(string message)	//If winners not found
	{
		throw message;
	}

	//Get the winner and weights
	weight.resize(testinput.size());
	for(unsigned int i=0; i<weight.size(); i++)
		weight[i] = somoutputlayer[position.x][position.y][position.z+i].weight;
	
	if(!outputpattern.empty())
		pattern = outputpattern[position.x][position.y];
	else
		pattern = -1;

	return;
}

MultiscaleNet::~MultiscaleNet(void)
{
}

void MultiscaleNet::SaveWeightsToFile(std::string filename)
{
	//Create a string to store the output of the weights
	std::ostringstream systemstring;
	systemstring << sizeofsom.x << '\t' << sizeofsom.y << "\r\n";
	for(unsigned int j=0; j<(unsigned int)somoutputlayer.size(); j++)
	{
		for(unsigned int i=0; i<(unsigned int)somoutputlayer[j].size(); i++)
		{
			for(unsigned int k=0; k<(unsigned int)somoutputlayer[j][i].size(); k++)
			{
				for(unsigned int l=0; l<(unsigned int)somoutputlayer[j][i][k].weight.size(); l++)
				{
					systemstring << somoutputlayer[j][i][k].weight[l] << '\t';
				}
			}
			systemstring << "\r\n";
		}
	}

	//Create a file stream to store the string
	FILE* file;
	file = fopen(filename.c_str(), "w");
	if(file == NULL)
		throw string("Error in opening file");
	else
	{
		if(fputs(systemstring.str().c_str(), file)==EOF)
			throw string("Error in writing file");
		fclose(file);
	}

	return;
}

_3dpoint MultiscaleNet::GetNetSize()
{
	return sizeofsom;
}

//Get the distance during training
vector<double> MultiscaleNet::GetDistances()
{
	if(trainingdistance.size()==0)
		throw string("SOM not trained");

	return trainingdistance;
}

//Save the net into file
void MultiscaleNet::SaveNetIntoFile(string filename)
{
	//Create a four dimension data for storing the weights
	vector<vector<vector<vector<double>>>> somdata;
	somdata.resize(sizeofsom.x);
	for(unsigned int i=0; i<sizeofsom.x; i++)
	{
		somdata[i].resize(sizeofsom.y);
		for(unsigned int j=0; j<sizeofsom.y; j++)
		{
			somdata[i][j].resize(sizeofsom.z);
			for(unsigned int k=0; k<sizeofsom.z; k++)
			{
				somdata[i][j][k] = somoutputlayer[i][j][k].weight;
			}
		}
	}

	//Save the data to file
	try
	{
		bf.SaveDataAsBinary(somdata, filename);
	}
	catch(string exceptionstring)
	{
		throw exceptionstring;
	}

	return;
}

//Load the net from a file
void MultiscaleNet::LoadNetFromFile(std::string filename)
{
	//Load the data from the file
	vector<vector<vector<vector<double>>>> somdata;
	try
	{
		bf.LoadBinaryData(somdata, filename);
	}
	catch(string exceptionstring)
	{
		throw exceptionstring;
	}

	//Load the data into the multiscale net
	sizeofsom.x = (unsigned int)somdata.size();
	sizeofsom.y = (unsigned int)somdata[0].size();
	sizeofsom.z = (unsigned int)somdata[0][0].size();
	datadimension = (unsigned int)somdata[0][0][0].size();

	somoutputlayer.resize(sizeofsom.x);
	for(unsigned int i=0; i<sizeofsom.x; i++)
	{
		somoutputlayer[i].resize(sizeofsom.y);
		if((unsigned int)somdata[i].size()!=sizeofsom.y)
			throw string("Data dimension not matched");		//Varify the loaded data
		for(unsigned int j=0; j<sizeofsom.y; j++)
		{
			somoutputlayer[i][j].resize(sizeofsom.z);
			if((unsigned int)somdata[i][j].size()!=sizeofsom.z)
				throw string("Data dimension not matched");		//Varify the loaded data
			for(unsigned int k=0; k<sizeofsom.z; k++)
			{
				if((unsigned int)somdata[i][j][k].size()!=datadimension)
					throw string("Data dimension not matched");		//Varify the loaded data
				somoutputlayer[i][j][k].weight = somdata[i][j][k];
				_3dpoint nodeposition = {i, j, k};
				somoutputlayer[i][j][k].position = nodeposition;
			}
		}
	}
	return;
}

//Get output layer patterns
vector<vector<unsigned int> > MultiscaleNet::GetOutputPattern()
{
	return outputpattern;
}
